Kongeørn

Analyser av kongeørnterritorier

Datasettet ligger på P:// som er en server på NINA.

#dat <- sf::st_read('P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/utbredelseKongeoern/Buff_10_All_terr_kongeørn.shp')


# Denne funker for rstudio server. Den over funker p lokal maskin.
dat <- sf::st_read('/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/utbredelseKongeoern/Buff_10_All_terr_kongeørn.shp')

## Reading layer `Buff_10_All_terr_kongeørn' from data source 
##   `/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/utbredelseKongeoern/Buff_10_All_terr_kongeørn.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1371 features and 3 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -60905 ymin: 6462156 xmax: 1102486 ymax: 7937378
## Projected CRS: WGS 84 / UTM zone 33N

head(dat)

## Simple feature collection with 6 features and 3 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 44492 ymin: 6590903 xmax: 79790 ymax: 6643075
## Projected CRS: WGS 84 / UTM zone 33N
##   Terr_total Visit_2010 Visited_20                       geometry
## 1  A-NAA-003          1          1 POLYGON ((74345 6633075, 73...
## 2  A-NAA-004          0          0 POLYGON ((79790 6624191, 79...
## 3  A-NAA-006          1          1 POLYGON ((72664 6617925, 72...
## 4  A-NAA-007          1          1 POLYGON ((64492 6602755, 64...
## 5  A-NAA-008          0          1 POLYGON ((72395 6607057, 71...
## 6  A-NAA-009          0          1 POLYGON ((74419 6600903, 73...

getwd()

## [1] "/data/Egenutvikling/41001581_egenutvikling_anders_kolstad/github/IBECA_NINA/R"

Datasettet inneholder alle kjente kongeørnterritorier i Norge (vi bryr oss ikke om streiffulger). Ikke alle territorier er okkupert til en hver tid. Det er to kolonner (Visit_2021 og Visit_20) som indikarer hvorvidt territoriene har blitt overvåket i hhv 2010-2014 og/eller 2015-2019. Det er med andre ord data mellom 2010 og 2019. La oss se hvor mange territorier det er:

length(unique(dat$Terr_total))

## [1] 1371

Er det noen duplikater?

nrow(dat)

## [1] 1371

Nei, det er det ikke.

Hvor mange territorier ble overvåket i siste periode?

table(dat$Visited_20)

## 
##   0   1 
## 838 533

Det er 533 territorier som ble overvåket i siste periode. Av disse var det kun 398 som var okkupert (pers. com. Jenny Mattison). Modelleringer tilseier at det skal være rundt 1027 okkuperte territorier i Norge under siste perioden (2015-2019) (pers. com. Jenny Mattison). Disse dataene finnes på fylkesnivå.

table(dat$Visit_2010, dat$Visited_20, deparse.level = 2)

##               dat$Visited_20
## dat$Visit_2010   0   1
##              0 635 198
##              1 203 335

Tabellen over sier oss at 635 territorier ikke har vært overvåket i noen av periodene (!?), og at 335 har vært overvåket i begge.

Ett territorie er en sirkel med radius 10 km plassert rundt sentrumspunktet for de reirene som er inkludert i det samme territoriet (det kan være flere, men bare ett er okkupert til en hvert tid).

La oss se hvordan territoriene fordeler seg i landet. Først, las oss se hvilke crs dataene våre har og deretter laste inn en outline for Norge.

crs(dat)

## CRS arguments:
##  +proj=utm +zone=33 +datum=WGS84 +units=m +no_defs

Outline of Norway

nor <- raster::getData('GADM', country='NOR', level=0)%>%
  sp::spTransform("+proj=utm +zone=33")
saveRDS(nor, '../data/norway_outline.RDS')

Lets also get a zoom box ready

myExt <- raster::extent(100000, 300000, 6900000, 7100000)

plot(nor, axes=T)
plot(dat[,c('Terr_total', 'geometry')], add=T)
plot(myExt, add=T, lwd=10)

Det er ingenting på Østlandet. Ellers så er landet godt dekt.

Zoomer inn i ruta…

datc <- sf::st_crop(dat, myExt)

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

plot(datc[,c('Terr_total', 'geometry')])

Det er som vi ser veldig mye overlapp mellom de teoretiske territoriene.

Fjellmasken

Det neste vi må gjøre er å ekskludere territorier som ikke har fjell. Her har vi flere muligheter. La oss først ta inn fjellmasken.

#fjell <- stars::read_stars("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/FINAL/Raw_data/Ecosystems/ecoMap_50m.tif",
#                           proxy=F)

#fjell <- raster::raster("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/FINAL/Raw_data/Ecosystems/ecoMap_50m.tif")
#fjell <- raster::raster("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/M/Fjell/Mountain ecosystem Norge 50m.tif")
fjell <- raster::raster("/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/M/Fjell/Mountain ecosystem Norge 50m.tif")

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on GRS80 ellipsoid in Proj4
## definition

class(fjell)

## [1] "RasterLayer"
## attr(,"package")
## [1] "raster"

Fila er ganske stor.

dim(fjell)

## [1] 32021 32776     1

32k * 32k gir litt over en milliard verdier (piksler)

plot(fjell)

Denne masken ser ut til å være piksler med verdi 202 fra dette kartlaget: P:\41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly_data_50m.tif

Fila er ganske tung å jobbe med så vi må prøve å redusere størrelsen

summary(fjell)

## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (0.01% of all cells)

##         Mountain_ecosystem_Norge_50m
## Min.                               0
## 1st Qu.                            0
## Median                             0
## 3rd Qu.                            0
## Max.                             202
## NA's                               0

Alle 0’er kan bli NA

# fjellc <-  raster::mask(fjell, fjell[fjell==0],maskvalue = 0, updatevalue = NA)

Opperasjonen over tar for lang tid (ike ferdig på 2 timer)

Finner centroiden til territoriene for å defeinere de til skog eller fjell

dat$centroid <- st_centroid(dat$geometry)

# crop again
datc <- sf::st_crop(dat, myExt)

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

datc$centroid <- st_centroid(datc$geometry)

nor_sf <- st_as_sf(nor)
norc <- sf::st_crop(nor_sf, myExt)

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

fjellc <- raster::crop(fjell, myExt)

fjellc[fjellc[]==0] <- NA
fjellc[fjellc[]>0] <- 1


par(mfrow=c(1,1))
plot(norc$geometry, axes=T)
plot(fjellc, add=T)
plot(datc$geometry, add=T)
plot(datc$centroid, add=T)

Om vi skal bruke centroidene til å kategoriese territoriene som skog- eller fjellterritorier så må nok fjellmasken være på den oppløsningen den er. Men om vi skal istede regne ut andel av fjellpiksler per territorie så kan vi redusere oppløsningen betraktelig uten å miste informasjon. Prøver først med centroidene:

datc2 <- as(datc$centroid, 'Spatial') # extract-funksjonen krever spatial points dataframe
datc2$class <- raster::extract(fjellc, datc2, fun=sum, df=F)

plot(norc$geometry, axes=T)
  plot(fjellc, add=T)
  plot(datc2[!is.na(datc2$class),], add=T, pch=21)
  plot(datc2[is.na(datc2$class),], add=T)

Her vises centroidene som ligger i fjell som fylte sirkler, og de som ligger utenfor som plusstegn. De fleste ligger utenfor. Dette er bare et lite utsnitt av Norge, Men det kan være representatbelt da kongeørn ofte vil ha trær for å bygge reir, men de kan fortsatt jakte på fjellet. Jeg prøver derfor også den andre metoden. Først reduserer jeg oppløsningen på fjellmasken fra 50x50m til 1000x1000 (factor = 20).

fjellc10 <- raster::aggregate(fjellc, fact = 20)
par(mfrow=c(1,2))
plot(fjellc); plot(fjellc10)

Noen av øyene ser ut til å ha blitt litt større. Kanskje vi kunne brukt fact = 10. La oss zoom inn litt til for å se.

myExt2 <- raster::extent(200000, 300000, 6900000, 7000000)
fjellcc <- crop(fjellc, myExt2)
fjellc10c <- crop(fjellc10, myExt2)
par(mfrow=c(1,2))
plot(fjellcc); plot(fjellc10c)

Dette kan være en konsekvens av at jeg satta 0-verdier til NA, slik at alle aggregerte celler blir 1 gitt at det er >=1 fjellcelle der.

fjellc[is.na(fjellc[])] <- 0
fjellc10v2 <- raster::aggregate(fjellc, fact = 20) # 5 sec. processing time
fjellccv2 <- crop(fjellc, myExt2)
fjellc10cv2 <- crop(fjellc10v2, myExt2)

# definerer nye fjellceller som de med fjellareal over 50%
fjellc10cv3 <- fjellc10cv2
fjellc10cv3[fjellc10cv3[]<0.5] <- 0
fjellc10cv3[fjellc10cv3[]>=0.5] <- 1

par(mfrow=c(2,2))
plot(fjellccv2, main="Original"); 
plot(fjellc10c, main = "With NA's instead of zeros"); 
plot(fjellc10cv2, main = "Aggregation with mean values"); 
plot(fjellc10cv3, main = "Aggregation with threshold")

Den siste ser grei ut. Plotter territoriene over masken for å sammenligne oppløsningen.

datcc <- st_crop(datc, myExt2)

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

plot(fjellc10cv3, main = "Aggregation with threshold")
  plot(datcc$geometry, add=T)

Ekstraktherer fraksjon fjellareal per polygon.

datcc$andelFjell <- raster::extract(fjellc10cv3, datcc, fun = mean) # 9 sec processing time

hist(datcc$andelFjell)

Territoriene har god spredning når det gjelder andel fjellareal

myCols <- sf.colors(5, alpha = 0.5)
datcc$colour <- myCols[1]
datcc$colour[datcc$andelFjell >=0.2 & datcc$andelFjell < 0.4] <- myCols[2]
datcc$colour[datcc$andelFjell >=0.4 & datcc$andelFjell < 0.6] <- myCols[3]
datcc$colour[datcc$andelFjell >=0.6 & datcc$andelFjell < 0.8] <- myCols[4]
datcc$colour[datcc$andelFjell >=0.8]                          <- myCols[5]

plot(fjellc10cv3, main = "Aggregation with threshold", legend=F)
  plot(datcc$geometry, add=T, col = datcc$colour)
  legend("bottomright",   
      legend = c("<.2", ".2-.4", ".4-.6", ".6-.8", ">.8"), 
      fill = myCols)

Dette differensierer territoriene etter hvor mye fjell de innehar. Gjennomsnittet i dette utsnittet er 0.4358142.

Sjekkliste

Sjekklisten er ikke utfylt enda.

Sjekkliste